!===========================================================================
! CAM interface to the UNIFIED CONVECTION SCHEME (UNICON)
!
! The USE_UNICON macro converts this module to a stub interface which allows
! CAM to be built without the unicon and unicon_utils modules.
!
!===========================================================================

module unicon_cam

  use shr_kind_mod    , only: r8 => shr_kind_r8, i4 => shr_kind_i4
  use spmd_utils      , only: masterproc
  use ppgrid          , only: pcols, pver, pverp, begchunk, endchunk
  use physconst       , only: rair, cpair, gravit, latvap, latice, zvir, mwdry
  use constituents    , only: pcnst, cnst_add, qmin, cnst_get_type_byind, cnst_get_ind, cnst_name
  use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_mode_num_idx, rad_cnst_get_mam_mmr_idx
  use physics_types   , only: physics_state, physics_ptend, physics_ptend_init
  use camsrfexch      , only: cam_in_t
  use physics_buffer  , only: pbuf_add_field, dtype_r8, dyn_time_lvls, pbuf_old_tim_idx, &
                              physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_set_field
  use phys_control    , only: phys_getopts
  use cam_history     , only: outfld, addfld, horiz_only
  use time_manager    , only: is_first_step
  use cam_abortutils  , only: endrun

#ifdef USE_UNICON
  use unicon          , only: unicon_init, compute_unicon
  use unicon_utils    , only: unicon_utils_init, positive_moisture, positive_tracer
#endif

  implicit none
  private
  save

  public &
    unicon_cam_readnl,      &
    unicon_cam_register,    &
    unicon_cam_init,        &
    unicon_implements_cnst, &
    unicon_init_cnst,       &
    unicon_out_t,           &
    unicon_cam_tend,        &
    unicon_cam_org_diags

  logical :: unicon_offline_dat_out   = .false.
  integer :: unicon_offline_dat_hfile = 2

  real(r8) xlv    ! Latent heat of vaporization
  real(r8) xlf    ! Latent heat of fusion
  real(r8) xls    ! Latent heat of sublimation
  real(r8) cp     ! Specific heat of dry air

  integer, parameter :: nseg = 1 ! Number of updraft segments [ # ]

  ! For advecting organization-related variables
  integer, parameter :: n_org = 5                  ! Number of constituents
  character(8), dimension(n_org), parameter :: &   ! Constituent names
    cnst_names = ['ORGawk  ','ORGthl  ','ORGqto  ','ORGuoo  ','ORGvoo  ']

  integer awk_cnst_ind, thl_cnst_ind, qt_cnst_ind, u_cnst_ind, v_cnst_ind

  ! Fields added to physics buffer by this module
  integer               &
    cushavg_idx       , &
    cuorg_idx         , &
    awk_PBL_idx       , &
    delta_thl_PBL_idx , &
    delta_qt_PBL_idx  , &
    delta_u_PBL_idx   , &
    delta_v_PBL_idx   , &
    delta_tr_PBL_idx  , &
    cu_cmfr_idx       , &
    cu_thlr_idx       , &
    cu_qtr_idx        , &
    cu_ur_idx         , &
    cu_vr_idx         , &
    cu_qlr_idx        , &
    cu_qir_idx        , &
    cu_trr_idx        , &
    cmfr_det_idx      , &
    qlr_det_idx       , &
    qir_det_idx

  ! Fields expected to be in the physics buffer
  integer ::              &
    sgh30_idx       = -1, &
    ast_idx         = -1, &
    tke_idx         = -1, &
    bprod_idx       = -1, &
    kpblh_idx       = -1, &
    pblh_idx        = -1, &
    went_idx        = -1, &
    cush_idx        = -1, &
    shfrc_idx       = -1, &
    icwmrsh_idx     = -1, &
    rprdsh_idx      = -1, &
    prec_sh_idx     = -1, &
    snow_sh_idx     = -1, &
    nevapr_shcu_idx = -1, &
    am_evp_st_idx   = -1, &   !  Evaporation area of stratiform precipitation [fraction]
    evprain_st_idx  = -1, &   !  Grid-mean evaporation rate of stratiform rain [kg/kg/s] >= 0.
    evpsnow_st_idx  = -1      !  Grid-mean evaporation rate of stratiform snow [kg/kg/s] >= 0.

  ! Constituent indices
  integer ixcldliq, ixcldice, ixnumliq, ixnumice

  ! Unicon output fields
  type unicon_out_t
    real(r8), allocatable, dimension(:,:) :: cmfmc ! Upward     convective mass flux at the interface [ kg / s / m2 ]
    real(r8), allocatable, dimension(:,:) :: slflx ! Net upward convective flux of liquid static energy [ J / s / m2 ]
    real(r8), allocatable, dimension(:,:) :: qtflx ! Net upward convective flux of total specific humidity [ kg / s / m2 ]
    real(r8), allocatable, dimension(:,:) :: rqc   ! Prod of suspended LWC+IWC by expel of excessive in-cumulus condensate [ kg / kg / s ] > 0
    real(r8), allocatable, dimension(:  ) :: rliq  ! Vertical integral of 'rqc' in flux unit [ m / s ]
    real(r8), allocatable, dimension(:  ) :: cnt   ! Cloud top  interface index ( ki = kpen )
    real(r8), allocatable, dimension(:  ) :: cnb   ! Cloud base interface index ( ki = krel-1 )
  contains
    procedure :: init  => unicon_out_init
    procedure :: clear => unicon_out_clear
    final unicon_out_final
  end type unicon_out_t

  ! Logical array to identify constituents that are mode number concentrations
  logical cnst_is_mam_num(PCNST)
  ! Logical array to identify constituents that are mode specie mass mixing ratios
  logical cnst_is_mam_mmr(PCNST)

contains

  subroutine unicon_out_init(this)

    class(unicon_out_t), intent(inout) :: this

    call this%clear()

    allocate(this%cmfmc(pcols,pver+1))
    allocate(this%slflx(pcols,pver+1))
    allocate(this%qtflx(pcols,pver+1))
    allocate(this%rqc  (pcols,pver  ))
    allocate(this%rliq (pcols       ))
    allocate(this%cnt  (pcols       ))
    allocate(this%cnb  (pcols       ))

  end subroutine unicon_out_init

  subroutine unicon_out_clear(this)

    class(unicon_out_t), intent(inout) :: this

    if (allocated(this%cmfmc)) deallocate(this%cmfmc)
    if (allocated(this%slflx)) deallocate(this%slflx)
    if (allocated(this%qtflx)) deallocate(this%qtflx)
    if (allocated(this%rqc  )) deallocate(this%rqc)
    if (allocated(this%rliq )) deallocate(this%rliq)
    if (allocated(this%cnt  )) deallocate(this%cnt)
    if (allocated(this%cnb  )) deallocate(this%cnb)

  end subroutine unicon_out_clear

  subroutine unicon_out_final(this)

    type(unicon_out_t), intent(inout) :: this

    call this%clear()

  end subroutine unicon_out_final

  subroutine unicon_cam_readnl(nlfile)

    use namelist_utils,  only: find_group_name
    use units,           only: getunit, freeunit
    use mpishorthand

    character(*), intent(in) :: nlfile

    integer unitn, ierr
    character(*), parameter :: subname = 'unicon_cam_readnl'

    namelist /unicon_nl/ unicon_offline_dat_out, unicon_offline_dat_hfile

    if (masterproc) then
      unitn = getunit()
      open(unitn, file=trim(nlfile), status='old')
      call find_group_name(unitn, 'unicon_nl', status=ierr)
      if (ierr == 0) then
        read(unitn, unicon_nl, iostat=ierr)
        if (ierr /= 0) then
          call endrun(subname // ':: ERROR reading namelist')
        end if
      end if
      close(unitn)
      call freeunit(unitn)
    end if

#ifdef SPMD
    ! Broadcast namelist variables
    call mpibcast (unicon_offline_dat_out,   1,  mpilog, 0, mpicom)
    call mpibcast (unicon_offline_dat_hfile, 1,  mpiint, 0, mpicom)
#endif

  end subroutine unicon_cam_readnl

  subroutine unicon_cam_register

    ! Jun.02.2012. Sungsu for advecting organization-related horizontal heterogeneity
    !              within PBL.
    !              For the time being, in order to save computation time, advection of aerosol perturbations
    !              are simply neglected.

#ifdef USE_UNICON
    call cnst_add(cnst_names(1), mwdry, cpair,    0._r8, awk_cnst_ind, &
       'Wake area within PBL associated with organization', readiv=.false., mixtype = 'dry')
    call cnst_add(cnst_names(2), mwdry, cpair,    0._r8, thl_cnst_ind, &
       'Perturbation of  thl associated with organization', readiv=.false., mixtype = 'dry')
    call cnst_add(cnst_names(3), mwdry, cpair,    0._r8,  qt_cnst_ind, &
       'Perturbation of  qt  associated with organization', readiv=.false., mixtype = 'dry')
    call cnst_add(cnst_names(4), mwdry, cpair,    0._r8,   u_cnst_ind, &
       'Perturbation of  u   associated with organization', readiv=.false., mixtype = 'dry')
    call cnst_add(cnst_names(5), mwdry, cpair,    0._r8,   v_cnst_ind, &
       'Perturbation of  v   associated with organization', readiv=.false., mixtype = 'dry')

    call pbuf_add_field('cushavg',       'global', dtype_r8, [pcols,dyn_time_lvls],            cushavg_idx)
    call pbuf_add_field('cuorg',         'global', dtype_r8, [pcols,dyn_time_lvls],            cuorg_idx)
    call pbuf_add_field('awk_PBL',       'global', dtype_r8, [pcols,dyn_time_lvls],            awk_PBL_idx)
    call pbuf_add_field('delta_thl_PBL', 'global', dtype_r8, [pcols,dyn_time_lvls],            delta_thl_PBL_idx)
    call pbuf_add_field('delta_qt_PBL',  'global', dtype_r8, [pcols,dyn_time_lvls],            delta_qt_PBL_idx)
    call pbuf_add_field('delta_u_PBL',   'global', dtype_r8, [pcols,dyn_time_lvls],            delta_u_PBL_idx)
    call pbuf_add_field('delta_v_PBL',   'global', dtype_r8, [pcols,dyn_time_lvls],            delta_v_PBL_idx)
    call pbuf_add_field('delta_tr_PBL',  'global', dtype_r8, [pcols,pcnst,dyn_time_lvls],      delta_tr_PBL_idx)
    call pbuf_add_field('cu_cmfr',       'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_cmfr_idx)
    call pbuf_add_field('cu_thlr',       'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_thlr_idx)
    call pbuf_add_field('cu_qtr',        'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_qtr_idx)
    call pbuf_add_field('cu_ur',         'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_ur_idx)
    call pbuf_add_field('cu_vr',         'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_vr_idx)
    call pbuf_add_field('cu_qlr',        'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_qlr_idx)
    call pbuf_add_field('cu_qir',        'global', dtype_r8, [pcols,pver,dyn_time_lvls],       cu_qir_idx)
    call pbuf_add_field('cu_trr',        'global', dtype_r8, [pcols,pver,pcnst,dyn_time_lvls], cu_trr_idx)
    call pbuf_add_field('cmfr_det',      'global', dtype_r8, [pcols,pver],                     cmfr_det_idx)
    call pbuf_add_field('qlr_det',       'global', dtype_r8, [pcols,pver],                     qlr_det_idx)
    call pbuf_add_field('qir_det',       'global', dtype_r8, [pcols,pver],                     qir_det_idx)
#endif

  end subroutine unicon_cam_register

  subroutine unicon_cam_init(pbuf2d)

    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    character(*), parameter :: sub='unicon_cam_init: '
    integer i, icnst, j, l, m, nmodes, nspec

    character(8 ) units
    character(30) varname
    character(60) surname
    character(2 ) numcha
    integer       msfc

#ifdef USE_UNICON
    xlv = latvap
    xlf = latice
    xls = xlv + xlf
    cp  = cpair

    call unicon_init(latvap, cpair, latice, zvir, rair, gravit)

    call cnst_get_ind('CLDLIQ', ixcldliq)
    call cnst_get_ind('CLDICE', ixcldice)
    call cnst_get_ind('NUMLIQ', ixnumliq)
    call cnst_get_ind('NUMICE', ixnumice)

    sgh30_idx       = pbuf_get_index('SGH30'      )
    ast_idx         = pbuf_get_index('AST'        )
    tke_idx         = pbuf_get_index('tke'        )
    bprod_idx       = pbuf_get_index('bprod'      )
    kpblh_idx       = pbuf_get_index('kpblh'      )
    pblh_idx        = pbuf_get_index('pblh'       )
    went_idx        = pbuf_get_index('went'       )
    cush_idx        = pbuf_get_index('cush'       )
    shfrc_idx       = pbuf_get_index('shfrc'      )
    icwmrsh_idx     = pbuf_get_index('ICWMRSH'    )
    rprdsh_idx      = pbuf_get_index('RPRDSH'     )
    prec_sh_idx     = pbuf_get_index('PREC_SH'    )
    snow_sh_idx     = pbuf_get_index('SNOW_SH'    )
    nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
    am_evp_st_idx   = pbuf_get_index('am_evp_st'  )
    evprain_st_idx  = pbuf_get_index('evprain_st' )
    evpsnow_st_idx  = pbuf_get_index('evpsnow_st' )

    ! physics buffer fields that need initializers -- these are only
    ! fields that have been added to pbuf by this module, or by the
    ! convection driver layer.
    if (is_first_step()) then
      if (cush_idx > 0) then
        call pbuf_set_field(pbuf2d, cush_idx,          1.e3_r8)
      else
        call endrun(sub // 'cush not in pbuf')
      end if
      if (cushavg_idx > 0) then
        call pbuf_set_field(pbuf2d, cushavg_idx,       1.e3_r8)
      else
        call endrun(sub//'cushavg not in pbuf')
      end if
      if (cuorg_idx > 0) then
        call pbuf_set_field(pbuf2d, cuorg_idx,         0.0_r8)
      else
        call endrun(sub//'cuorg not in pbuf')
      end if
      if (awk_pbl_idx > 0) then
        call pbuf_set_field(pbuf2d, awk_pbl_idx,       0.0_r8)
      else
        call endrun(sub//'awk_PBL not in pbuf')
      end if
      if (delta_thl_PBL_idx > 0) then
        call pbuf_set_field(pbuf2d, delta_thl_PBL_idx, 0.0_r8)
      else
        call endrun(sub//'delta_thl_PBL not in pbuf')
      end if
      if (delta_qt_PBL_idx > 0) then
        call pbuf_set_field(pbuf2d, delta_qt_PBL_idx,  0.0_r8)
      else
        call endrun(sub//'delta_qt_PBL not in pbuf')
      end if
      if (delta_u_PBL_idx > 0) then
        call pbuf_set_field(pbuf2d, delta_u_PBL_idx,   0.0_r8)
      else
        call endrun(sub//'delta_u_PBL not in pbuf')
      end if
      if (delta_v_PBL_idx > 0) then
        call pbuf_set_field(pbuf2d, delta_v_PBL_idx,   0.0_r8)
      else
        call endrun(sub//'delta_v_PBL not in pbuf')
      end if
      if (delta_tr_PBL_idx > 0) then
        call pbuf_set_field(pbuf2d, delta_tr_PBL_idx,  0.0_r8)
      else
        call endrun(sub//'delta_tr_PBL not in pbuf')
      end if
      if (cu_cmfr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_cmfr_idx,       0.0_r8)
      else
        call endrun(sub//'cu_cmfr not in pbuf')
      end if
      if (cu_thlr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_thlr_idx,       0.0_r8)
      else
        call endrun(sub//'cu_thlr not in pbuf')
      end if
      if (cu_qtr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_qtr_idx,        0.0_r8)
      else
        call endrun(sub//'cu_qtr not in pbuf')
      end if
      if (cu_ur_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_ur_idx,         0.0_r8)
      else
        call endrun(sub//'cu_ur not in pbuf')
      end if
      if (cu_vr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_vr_idx,         0.0_r8)
      else
        call endrun(sub//'cu_vr not in pbuf')
      end if
      if (cu_qlr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_qlr_idx,        0.0_r8)
      else
        call endrun(sub//'cu_qlr not in pbuf')
      end if
      if (cu_qir_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_qir_idx,        0.0_r8)
      else
        call endrun(sub//'cu_qir not in pbuf')
      end if
      if (cu_trr_idx > 0) then
        call pbuf_set_field(pbuf2d, cu_trr_idx,        0.0_r8)
      else
        call endrun(sub//'cu_trr not in pbuf')
      end if
      if (cmfr_det_idx > 0) then
        call pbuf_set_field(pbuf2d, cmfr_det_idx,      0.0_r8)
      else
        call endrun(sub//'cmfr_det not in pbuf')
      end if
      if (qlr_det_idx > 0) then
        call pbuf_set_field(pbuf2d, qlr_det_idx,      0.0_r8)
      else
        call endrun(sub//'qlr_det not in pbuf')
      end if
      if (qir_det_idx > 0) then
        call pbuf_set_field(pbuf2d, qir_det_idx,      0.0_r8)
      else
        call endrun(sub//'qir_det not in pbuf')
      end if
    end if

    ! Set arrays to identify the modal aerosol constituents
    cnst_is_mam_num = .false.
    cnst_is_mam_mmr = .false.
    call rad_cnst_get_info(0, nmodes=nmodes)
    do i = 1, nmodes
      call rad_cnst_get_mode_num_idx(i, icnst)
      cnst_is_mam_num(icnst) = .true.
      call rad_cnst_get_info(0, i, nspec=nspec)
      do j = 1, nspec
        call rad_cnst_get_mam_mmr_idx(i, j, icnst)
        cnst_is_mam_mmr(icnst) = .true.
      end do
    end do

    ! ------------------------- !
    ! Internal Output Variables !
    ! ------------------------- !

    ! Sungsu for advection of convective organization

    call addfld('ORGawk', [ 'lev' ],  'A', 'no'   , 'Convective Organization - Wake Area' )
    call addfld('ORGthl', [ 'lev' ],  'A', 'K'    , 'Convective Organization - Perturbation of thl in the non-wake area')
    call addfld('ORGqto', [ 'lev' ],  'A', 'kg/kg', 'Convective Organization - Perturbation of qt  in the non-wake area')
    call addfld('ORGuoo', [ 'lev' ],  'A', 'm/s'  , 'Convective Organization - Perturbation of u  in the non-wake area' )
    call addfld('ORGvoo', [ 'lev' ],  'A', 'm/s'  , 'Convective Organization - Perturbation of v  in the non-wake area' )

    ! From the main unified convection scheme

    call addfld('flxrain_SP'  , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective net rain flux')
    call addfld('flxsnow_SP'  , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective net snow flux')
    call addfld('cmf_SP'      , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective net mass flux')
    call addfld('qtflx_SP'    , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective net qt flux')
    call addfld('slflx_SP'    , [ 'ilev' ], 'A', 'J/m2/s' , 'Convective net sl flux')
    call addfld('uflx_SP'     , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective net u flux')
    call addfld('vflx_SP'     , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective net v flux')
    call addfld('cmf_u_SP'    , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective updraft mass flux')
    call addfld('qtflx_u_SP'  , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective updraft qt flux')
    call addfld('slflx_u_SP'  , [ 'ilev' ], 'A', 'J/m2/s' , 'Convective updraft sl flux')
    call addfld('uflx_u_SP'   , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective updraft u flux')
    call addfld('vflx_u_SP'   , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective updraft v flux')
    call addfld('cmf_d_SP'    , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective downdraft mass flux')
    call addfld('qtflx_d_SP'  , [ 'ilev' ], 'A', 'kg/m2/s', 'Convective downdraft qt flux')
    call addfld('slflx_d_SP'  , [ 'ilev' ], 'A', 'J/m2/s' , 'Convective downdraft sl flux')
    call addfld('uflx_d_SP'   , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective downdraft u flux')
    call addfld('vflx_d_SP'   , [ 'ilev' ], 'A', 'kg/m/s2', 'Convective downdraft v flux')
    call addfld('qtten_u_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qt by updraft')
    call addfld('slten_u_SP'  , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency sl by updraft')
    call addfld('uten_u_SP'   , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency u  by updraft')
    call addfld('vten_u_SP'   , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency v  by updraft')
    call addfld('sten_u_SP'   , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency s  by updraft')
    call addfld('qvten_u_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qv by updraft')
    call addfld('qlten_u_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by updraft')
    call addfld('qiten_u_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by updraft')
    call addfld('nlten_u_SP'  , [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency nl by updraft')
    call addfld('niten_u_SP'  , [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency ni by updraft')
    call addfld('qtten_d_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qt by downdraft')
    call addfld('slten_d_SP'  , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency sl by downdraft')
    call addfld('uten_d_SP'   , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency u  by downdraft')
    call addfld('vten_d_SP'   , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency v  by downdraft')
    call addfld('sten_d_SP'   , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency s  by downdraft')
    call addfld('qvten_d_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qv by downdraft')
    call addfld('qlten_d_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by downdraft')
    call addfld('qiten_d_SP'  , [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by downdraft')
    call addfld('nlten_d_SP'  , [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency nl by downdraft')
    call addfld('niten_d_SP'  , [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency ni by downdraft')
    call addfld('qtten_evp_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qt by evap of precip within environment')
    call addfld('slten_evp_SP', [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency sl by evap of precip within environment')
    call addfld('uten_evp_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency u  by evap of precip within environment')
    call addfld('vten_evp_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency v  by evap of precip within environment')
    call addfld('sten_evp_SP' , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency s  by evap of precip within environment')
    call addfld('qvten_evp_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qv by evap of precip within environment')
    call addfld('qlten_evp_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by evap of precip within environment')
    call addfld('qiten_evp_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by evap of precip within environment')
    call addfld('nlten_evp_SP', [ 'lev'  ], 'A', '#/kg/s' , 'Convective tendency nl by evap of precip within environment')
    call addfld('niten_evp_SP', [ 'lev'  ], 'A', '#/kg/s' , 'Convective tendency ni by evap of precip within environment')
    call addfld('qtten_dis_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qt by dissipative heating')
    call addfld('slten_dis_SP', [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency sl by dissipative heating')
    call addfld('uten_dis_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency u  by dissipative heating')
    call addfld('vten_dis_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency v  by dissipative heating')
    call addfld('sten_dis_SP' , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency s  by dissipative heating')
    call addfld('qvten_dis_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qv by dissipative heating')
    call addfld('qlten_dis_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by dissipative heating')
    call addfld('qiten_dis_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by dissipative heating')
    call addfld('nlten_dis_SP', [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency nl by dissipative heating')
    call addfld('niten_dis_SP', [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency ni by dissipative heating')
    call addfld('qtten_pos_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qt by positive tracer constraints')
    call addfld('slten_pos_SP', [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency sl by positive tracer constraints')
    call addfld('uten_pos_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency u  by positive tracer constraints')
    call addfld('vten_pos_SP' , [ 'lev'  ], 'A', 'm/s/s'  , 'Convective tendency v  by positive tracer constraints')
    call addfld('sten_pos_SP' , [ 'lev'  ], 'A', 'J/kg/s' , 'Convective tendency s  by positive tracer constraints')
    call addfld('qvten_pos_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qv by positive tracer constraints')
    call addfld('qlten_pos_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by positive tracer constraints')
    call addfld('qiten_pos_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by positive tracer constraints')
    call addfld('nlten_pos_SP', [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency nl by positive tracer constraints')
    call addfld('niten_pos_SP', [ 'lev'  ], 'A', '1/kg/s' , 'Convective tendency ni by positive tracer constraints')
    call addfld('qlten_sub_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by compensating subsidence')
    call addfld('qiten_sub_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by compensating subsidence')
    call addfld('qlten_det_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency ql by detrainment')
    call addfld('qiten_det_SP', [ 'lev'  ], 'A', 'kg/kg/s', 'Convective tendency qi by detrainment')
    call addfld('CMFR_DET'    , [ 'lev'  ], 'A', 'kg/m2/s', 'Detrained convective mass flux'       )
    call addfld('QLR_DET'     , [ 'lev'  ], 'A', 'kg/kg'  , 'Detrained convective LWC'             )
    call addfld('QIR_DET'     , [ 'lev'  ], 'A', 'kg/kg'  , 'Detrained convective IWC'             )
    call addfld('cush_SP'     , horiz_only, 'A', 'm'      , 'Cumulus top height from UNICON')
    call addfld('cushavg_SP'  , horiz_only, 'A', 'm'      , 'Mean cumulus top height from UNICON')
    call addfld('cuorg_SP'    , horiz_only, 'A', '1'      , 'Cumulus organization parameter from UNICON')
    call addfld('Radius_SP'   , horiz_only, 'A', 'm'      , 'Cumulus plume radius at surface from UNICON')
    call addfld('sgh30_SP'    , horiz_only, 'A', 'm'      , 'Standard deviation of subgrid topography at 30 sec')
    call addfld('kw_SP'       , horiz_only, 'A', 'no'     , 'Internally computed kw from SPARKCONV' )
    call addfld('sigma_w_SP'  , horiz_only, 'A', 'm/s'    , 'Standard deviation of updraft w at surface from UNICON')
    call addfld('sigma_thl_SP', horiz_only, 'A', 'K'      , 'Standard deviation of updraft thl at surface from UNICON')
    call addfld('sigma_qt_SP' , horiz_only, 'A', 'g/kg'   , 'Standard deviation of updraft qt at surface from UNICON')
    call addfld('sigma_u_SP'  , horiz_only, 'A', 'm/s'    , 'Standard deviation of updraft u at surface from UNICON')
    call addfld('sigma_v_SP'  , horiz_only, 'A', 'm/s'    , 'Standard deviation of updraft v at surface from UNICON')
    call addfld('w_org_SP'    , horiz_only, 'A', 'm2/s2'  , 'Organization-generated additional w at surface from UNICON')
    call addfld('thl_org_SP'  , horiz_only, 'A', 'K'      , 'Organization-generated additional thl at surface from UNICON')
    call addfld('qt_org_SP'   , horiz_only, 'A', 'g/kg'   , 'Organization-generated additional qt at surface from UNICON')
    call addfld('u_org_SP'    , horiz_only, 'A', 'm/s'    , 'Organization-generated additional u at surface from UNICON')
    call addfld('v_org_SP'    , horiz_only, 'A', 'm/s'    , 'Organization-generated additional v at surface from UNICON')
    call addfld('tkes_SP'     , horiz_only, 'A', 'm2/s2'  , 'TKE at surface within UNICON')
    call addfld('went_SP'     , horiz_only, 'A', 'm/s'    , 'Entrainment rate at the PBL top interface from UW PBL')
    call addfld('went_eff_SP' , horiz_only, 'A', 'm/s'    , 'Effective entrainment rate at the PBL top interface in UNICON')
    call addfld('am_u_SP'     , [ 'lev'  ], 'A', '1'      , 'Convective updraft fractional area at mid-layer')
    call addfld('am_d_SP'     , [ 'lev'  ], 'A', '1'      , 'Convective downdraft fractional area at mid-layer')
    call addfld('qlm_u_SP'    , [ 'lev'  ], 'A', 'kg/kg'  , 'Area-weighted in-cumulus LWC condensate of updraft at mid-layer')
    call addfld('qlm_d_SP'    , [ 'lev'  ], 'A', 'kg/kg'  , 'Area-weighted in-cumulus IWC condensate of downdraft at mid-layer')
    call addfld('qim_u_SP'    , [ 'lev'  ], 'A', 'kg/kg'  , 'Area-weighted in-cumulus LWC condensate of updraft at mid-layer')
    call addfld('qim_d_SP'    , [ 'lev'  ], 'A', 'kg/kg'  , 'Area-weighted in-cumulus IWC condensate of downdraft at mid-layer')
    call addfld('thl_u_SP'    , [ 'ilev' ], 'A', 'K'      , 'Mass-flux weighted updraft mean thl')
    call addfld('qt_u_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted updraft mean qt')
    call addfld('u_u_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted updraft mean u')
    call addfld('v_u_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted updraft mean v')
    call addfld('w_u_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted updraft mean w')
    call addfld('ql_u_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted updraft mean ql')
    call addfld('qi_u_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted updraft mean qi')
    call addfld('wa_u_SP'     , [ 'ilev' ], 'A', 'm/s'    , 'Area weighted updraft mean w')
    call addfld('qla_u_SP'    , [ 'ilev' ], 'A', 'kg/kg'  , 'Area weighted updraft mean ql')
    call addfld('qia_u_SP'    , [ 'ilev' ], 'A', 'kg/kg'  , 'Area weighted updraft mean qi')
    call addfld('a_u_SP'      , [ 'ilev' ], 'A', '1'      , 'Convective updraft fractional area')
    call addfld('rad_u_SP'    , [ 'ilev' ], 'A', 'm'      , 'Number-weighted effective radius of updraft plumes')
    call addfld('num_u_SP'    , [ 'ilev' ], 'A', '1/m^2'  , 'Number concentration of updraft plumes')
    call addfld('gamw_u_SP'   , [ 'ilev' ], 'A', 'ratio'  , 'The ratio of w_u to wa_u')
    call addfld('nl_u_SP'     , [ 'ilev' ], 'A', '1/kg'   , 'Mass-flux weighted updraft mean nl')
    call addfld('ni_u_SP'     , [ 'ilev' ], 'A', '1/kg'   , 'Mass-flux weighted updraft mean ni')
    call addfld('thva_u_SP'   , [ 'ilev' ], 'A', 'K'      , 'Area weighted updraft mean thv')
    call addfld('thl_d_SP'    , [ 'ilev' ], 'A', 'K'      , 'Mass-flux weighted downdraft mean thl')
    call addfld('qt_d_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted downdraft mean qt')
    call addfld('u_d_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted downdraft mean u')
    call addfld('v_d_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted downdraft mean v')
    call addfld('w_d_SP'      , [ 'ilev' ], 'A', 'm/s'    , 'Mass-flux weighted downdraft mean w')
    call addfld('ql_d_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted downdraft mean ql')
    call addfld('qi_d_SP'     , [ 'ilev' ], 'A', 'kg/kg'  , 'Mass-flux weighted downdraft mean qi')
    call addfld('wa_d_SP'     , [ 'ilev' ], 'A', 'm/s'    , 'Area weighted downdraft mean w')
    call addfld('qla_d_SP'    , [ 'ilev' ], 'A', 'kg/kg'  , 'Area weighted downdraft mean ql')
    call addfld('qia_d_SP'    , [ 'ilev' ], 'A', 'kg/kg'  , 'Area weighted downdraft mean qi')
    call addfld('a_d_SP'      , [ 'ilev' ], 'A', '1'      , 'Convective downdraft fractional area')
    call addfld('nl_d_SP'     , [ 'ilev' ], 'A', '1/kg'   , 'Mass-flux weighted downdraft mean nl')
    call addfld('ni_d_SP'     , [ 'ilev' ], 'A', '1/kg'   , 'Mass-flux weighted downdraft mean ni')
    call addfld('thv_b_SP'    , [ 'ilev' ], 'A', 'K'      , 'thv_b : Environmental buoyancy at the base interface')
    call addfld('thv_t_SP'    , [ 'ilev' ], 'A', 'K'      , 'thv_t : Environmental buoyancy at the top interface')
    call addfld('thv_mt_SP'   , [ 'ilev' ], 'A', 'K'      , 'thv_mt : Environmental buoyancy at the top interface of lower layer')
    call addfld('thv_min_SP'  , [ 'ilev' ], 'A', 'K'      , 'thv_min : Minimum environmental buoyancy for downdraft sorting')
    call addfld('cu_cmfr_SP'  , [ 'lev'  ], 'A', 'kg/m2/s', 'Mass flux of mixing environmental airs')
    call addfld('cu_thlr_SP'  , [ 'lev'  ], 'A', 'K'      , 'thl       of mixing environmental airs')
    call addfld('cu_qtr_SP'   , [ 'lev'  ], 'A', 'kg/kg'  , 'qt        of mixing environmental airs')
    call addfld('cu_qlr_SP'   , [ 'lev'  ], 'A', 'kg/kg'  , 'ql        of mixing environmental airs')
    call addfld('cu_qir_SP'   , [ 'lev'  ], 'A', 'kg/kg'  , 'qi        of mixing environmental airs')
    call addfld('cu_ur_SP'    , [ 'lev'  ], 'A', 'm/s'    , 'u         of mixing environmental airs')
    call addfld('cu_vr_SP'    , [ 'lev'  ], 'A', 'm/s'    , 'v         of mixing environmental airs')
    call addfld('cu_thvr_SP'  , [ 'lev'  ], 'A', 'K'      , 'thv       of mixing environmental airs')
    call addfld('cu_rhr_SP'   , [ 'lev'  ], 'A', '1'      , 'rh        of mixing environmental airs')
    call addfld('cu_nlr_SP'   , [ 'lev'  ], 'A', '1/kg'   , 'nl        of mixing environmental airs')
    call addfld('cu_nir_SP'   , [ 'lev'  ], 'A', '1/kg'   , 'ni        of mixing environmental airs')
    call addfld('a_p_SP'      , [ 'ilev' ], 'A', '1'      , 'Convective precipitation area')
    call addfld('am_evp_SP'   , [ 'lev'  ], 'A', '1'      , 'Convective evaporation area')
    call addfld('am_pu_SP'    , [ 'lev'  ], 'A', '1'      , 'Overlapping area between conv precipitation and sat updraft area')
    call addfld('x_um_SP'     , [ 'lev'  ], 'A', 'm'      , 'Zonal displacement of the updraft area from the surface')
    call addfld('y_um_SP'     , [ 'lev'  ], 'A', 'm'      , 'Meridional displacement of the updraft area from the surface')
    call addfld('x_p_SP'      , [ 'ilev' ], 'A', 'm'      , 'Zonal displacement of the precipitation area from the surface')
    call addfld('y_p_SP'      , [ 'ilev' ], 'A', 'm'      , 'Meridional displacement of the precipitation area from the surface')

    do l = 1, pcnst
      if (cnst_is_mam_num(l) .or. cnst_is_mam_mmr(l)) then
         units = '1/kg/s'
         if (cnst_is_mam_mmr(l)) units = 'kg/kg/s'

         varname = trim(cnst_name(l))//'_u_SP'
         surname = trim(cnst_name(l))//' tendency by convective updraft from UNICON'
         call addfld(trim(varname), [ 'lev' ], 'A', units, trim(surname))

         varname = trim(cnst_name(l))//'_d_SP'
         surname = trim(cnst_name(l))//' tendency by convective downdraft from UNICON'
         call addfld(trim(varname), [ 'lev' ], 'A', units, trim(surname))

         varname = trim(cnst_name(l))//'_evp_SP'
         surname = trim(cnst_name(l))//' tendency by evap. of precip in env. from UNICON'
         call addfld(trim(varname), [ 'lev' ], 'A', units, trim(surname))

         varname = trim(cnst_name(l))//'_dis_SP'
         surname = trim(cnst_name(l))//' tendency by dissipative heating from UNICON'
         call addfld(trim(varname), [ 'lev' ], 'A', units, trim(surname))

         varname = trim(cnst_name(l))//'_pos_SP'
         surname = trim(cnst_name(l))//' tendency by positive moisture from UNICON'
         call addfld(trim(varname), [ 'lev' ], 'A', units, trim(surname))
      end if
    end do

    ! Nov.15.2012. Below output corresponding to individual updraft segment is designed to write out individual
    !              segment values for writing UNICON_II paper.
    do msfc = 1, nseg
      write(numcha, '(i2.2)') msfc
      ! The properties of individual updraft segment
      call addfld('thl_u'     //numcha//'_SP', ['ilev'], 'A', 'K'       , numcha//' updraft segment : updraft thl'            )
      call addfld('qt_u'      //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , numcha//' updraft segment : updraft qt'             )
      call addfld('u_u'       //numcha//'_SP', ['ilev'], 'A', 'm/s'     , numcha//' updraft segment : updraft u'              )
      call addfld('v_u'       //numcha//'_SP', ['ilev'], 'A', 'm/s'     , numcha//' updraft segment : updraft v'              )
      call addfld('w_u'       //numcha//'_SP', ['ilev'], 'A', 'm/s'     , numcha//' updraft segment : updraft w'              )
      call addfld('ql_u'      //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , numcha//' updraft segment : updraft ql'             )
      call addfld('qi_u'      //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , numcha//' updraft segment : updraft qi'             )
      call addfld('cmf_u'     //numcha//'_SP', ['ilev'], 'A', 'kg/s/m^2', numcha//' updraft segment : updraft cmf'            )
      call addfld('a_u'       //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft fractional area')
      call addfld('num_u'     //numcha//'_SP', ['ilev'], 'A', '1/m^2'   , numcha//' updraft segment : updraft number density' )
      call addfld('rad_u'     //numcha//'_SP', ['ilev'], 'A', 'm'       , numcha//' updraft segment : updraft plume radius'   )
      call addfld('nl_u'      //numcha//'_SP', ['ilev'], 'A', '1/kg'    , numcha//' updraft segment : updraft nl'             )
      call addfld('ni_u'      //numcha//'_SP', ['ilev'], 'A', '1/kg'    , numcha//' updraft segment : updraft ni'             )
      call addfld('eps0_u'    //numcha//'_SP', ['ilev'], 'A', '1/Pa'    , numcha//' updraft segment : updraft eps0'    )
      call addfld('eps_u'     //numcha//'_SP', ['ilev'], 'A', '1/Pa'    , numcha//' updraft segment : updraft eps'     )
      call addfld('del_u'     //numcha//'_SP', ['ilev'], 'A', '1/Pa'    , numcha//' updraft segment : updraft del'     )
      call addfld('eeps_u'    //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft eeps'    )
      call addfld('ddel_u'    //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft ddel'    )
      call addfld('xc_u'      //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft xc'      )
      call addfld('xs_u'      //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft xs'      )
      call addfld('xemin_u'   //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft xemin'   )
      call addfld('xemax_u'   //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft xemax'   )
      call addfld('cridis_u'  //numcha//'_SP', ['ilev'], 'A', 'm'       , numcha//' updraft segment : updraft cridis'  )
      call addfld('thvcuenv_u'//numcha//'_SP', ['ilev'], 'A', 'K'       , numcha//' updraft segment : updraft thvcuenv')
      call addfld('thvegenv_u'//numcha//'_SP', ['ilev'], 'A', 'K'       , numcha//' updraft segment : updraft thvegenv')
      call addfld('thvxsenv_u'//numcha//'_SP', ['ilev'], 'A', 'K'       , numcha//' updraft segment : updraft thvxsenv')
      call addfld('fmix_u'    //numcha//'_SP', ['ilev'], 'A', '1'       , numcha//' updraft segment : updraft fmix'    )
      call addfld('cmfumix_u' //numcha//'_SP', ['ilev'], 'A', 'kg/s/m^2', numcha//' updraft segment : updraft cmfumix' )
      call addfld('ptop'      //numcha//'_SP', horiz_only, 'A', 'Pa', numcha//' updraft segment : updraft top pressure')
      call addfld('ztop'      //numcha//'_SP', horiz_only, 'A', 'm' , numcha//' updraft segment : updraft top height'  )

      ! The properties of mass flux weighted ( or area-weighted or net=sum ) downdraft properties for individual updraft segment

      call addfld('thl_d' //numcha//'_SP', ['ilev'], 'A', 'K'       , 'Mass-flux weighted  mean downdraft thl for '// numcha//' updraft segment')
      call addfld('qt_d'  //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , 'Mass-flux weighted  mean downdraft  qt for '// numcha//' updraft segment')
      call addfld('u_d'   //numcha//'_SP', ['ilev'], 'A', 'm/s'     , 'Mass-flux weighted  mean downdraft   u for '// numcha//' updraft segment')
      call addfld('v_d'   //numcha//'_SP', ['ilev'], 'A', 'm/s'     , 'Mass-flux weighted  mean downdraft   v for '// numcha//' updraft segment')
      call addfld('w_d'   //numcha//'_SP', ['ilev'], 'A', 'm/s'     , 'Mass-flux weighted  mean downdraft   w for '// numcha//' updraft segment')
      call addfld('ql_d'  //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , 'Mass-flux weighted  mean downdraft  ql for '// numcha//' updraft segment')
      call addfld('qi_d'  //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , 'Mass-flux weighted  mean downdraft  qi for '// numcha//' updraft segment')
      call addfld('wa_d'  //numcha//'_SP', ['ilev'], 'A', 'm/s'     , 'Area      weighted  mean downdraft   w for '// numcha//' updraft segment')
      call addfld('qla_d' //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , 'Area      weighted  mean downdraft  ql for '// numcha//' updraft segment')
      call addfld('qia_d' //numcha//'_SP', ['ilev'], 'A', 'kg/kg'   , 'Area      weighted  mean downdraft  qi for '// numcha//' updraft segment')
      call addfld('cmf_d' //numcha//'_SP', ['ilev'], 'A', 'kg/s/m^2', 'Net                      downdraft cmf for '// numcha//' updraft segment')
      call addfld('a_d'   //numcha//'_SP', ['ilev'], 'A', 'fraction', 'Net                      downdraft   a for '// numcha//' updraft segment')
      call addfld('nl_d'  //numcha//'_SP', ['ilev'], 'A', '#/kg'    , 'Mass-flux weighted  mean downdraft  nl for '// numcha//' updraft segment')
      call addfld('ni_d'  //numcha//'_SP', ['ilev'], 'A', '#/kg'    , 'Mass-flux weighted  mean downdraft  ni for '// numcha//' updraft segment')
    end do

    ! Nov.16.2012. Additional detailed diagnostic output
    call addfld('thl_orgfce_SP'  , horiz_only, 'A', 'K/s'     , 'Total organization forcing generating thl difference between non-wake and grid-mean areas')
    call addfld('qt_orgfce_SP'   , horiz_only, 'A', 'kg/kg/s' , 'Total organization forcing generating qt  difference between non-wake and grid-mean areas')
    call addfld('u_orgfce_SP'    , horiz_only, 'A', 'm/s/s'   , 'Total organization forcing generating u difference between non-wake and grid-mean areas')
    call addfld('v_orgfce_SP'    , horiz_only, 'A', 'm/s/s'   , 'Total organization forcing generating v difference between non-wake and grid-mean areas')
    call addfld('awk_orgfce_SP'  , horiz_only, 'A', '1/s'     , 'Total organization forcing generating wake area')
    call addfld('thl_orgfce_f_SP', horiz_only, 'A', 'K/s'     , 'PBL top flux-related forcing generating thl difference between non-wake and grid-mean areas')
    call addfld('qt_orgfce_f_SP' , horiz_only, 'A', 'kg/kg/s' , 'PBL top flux-related forcing generating qt difference between non-wake and grid-mean areas')
    call addfld('u_orgfce_f_SP'  , horiz_only, 'A', 'm/s/s'   , 'PBL top flux-related forcing generating u difference between non-wake and grid-mean areas')
    call addfld('v_orgfce_f_SP'  , horiz_only, 'A', 'm/s/s'   , 'PBL top flux-related forcing generating v difference between non-wake and grid-mean areas')
    call addfld('awk_orgfce_f_SP', horiz_only, 'A', '1/s'     , 'PBL top flux-related forcing generating wake area')
    call addfld('thl_orgfce_u_SP', horiz_only, 'A', 'K/s'     , 'Up-and-Down diabatic forcing generating thl difference between non-wake and grid-mean areas')
    call addfld('qt_orgfce_u_SP' , horiz_only, 'A', 'kg/kg/s' , 'Up-and-Down diabatic forcing generating qt difference between non-wake and grid-mean areas')
    call addfld('u_orgfce_u_SP'  , horiz_only, 'A', 'm/s/s'   , 'Up-and-Down diabatic forcing generating u difference between non-wake and grid-mean areas')
    call addfld('v_orgfce_u_SP'  , horiz_only, 'A', 'm/s/s'   , 'Up-and-Down diabatic forcing generating v difference between non-wake and grid-mean areas')
    call addfld('awk_orgfce_m_SP', horiz_only, 'A', '1/s'     , 'Lateral-Mixing forcing for wake area')
    call addfld('thl_orgfce_e_SP', horiz_only, 'A', 'K/s'     , 'Environment diabatic forcing generating thl difference between non-wake and grid-mean areas')
    call addfld('qt_orgfce_e_SP' , horiz_only, 'A', 'kg/kg/s' , 'Environment diabatic forcing generating qt difference between non-wake and grid-mean areas')
    call addfld('u_orgfce_e_SP'  , horiz_only, 'A', 'm/s/s'   , 'Environment diabatic forcing generating u difference between non-wake and grid-mean areas')
    call addfld('v_orgfce_e_SP'  , horiz_only, 'A', 'm/s/s'   , 'Environment diabatic forcing generating v difference between non-wake and grid-mean areas')
    call addfld('cmf_d_orgh_SP'  , horiz_only, 'A', 'kg/m^2/s', 'Organization-inducing downdraft mass flux at the PBL top interface')
    call addfld('taui_thl_SP'    , horiz_only, 'A', '1/s'     , 'Inverse of damping time scale of the difference between off-wake and grid-mean thl')
    call addfld('taui_qt_SP'     , horiz_only, 'A', '1/s'     , 'Inverse of damping time scale of the difference between off-wake and grid-mean qt' )
    call addfld('taui_u_SP'      , horiz_only, 'A', '1/s'     , 'Inverse of damping time scale of the difference between off-wake and grid-mean u'  )
    call addfld('taui_v_SP'      , horiz_only, 'A', '1/s'     , 'Inverse of damping time scale of the difference between off-wake and grid-mean v'  )
    call addfld('taui_awk_SP'    , horiz_only, 'A', '1/s'     , 'Inverse of damping time scale of the wake area')
    call addfld('del_org_SP'     , horiz_only, 'A', '1/s'     , 'Detrainment rate of the cold pool from UNICON')
    call addfld('del0_org_SP'    , horiz_only, 'A', '1/s'     , 'Effective detrainment rate of the cold pool from UNICON')
#endif

  end subroutine unicon_cam_init

  logical function unicon_implements_cnst(name) result(res)

    character(*), intent(in) :: name  ! constituent name

    integer m

    res = .false.

#ifdef USE_UNICON
    do m = 1, n_org
      if (name == cnst_names(m)) then
        res = .true.
        return
      end if
    end do
#endif

  end function unicon_implements_cnst

  subroutine unicon_init_cnst(name, latvals, lonvals, mask, q)

    character(*), intent(in ) :: name       ! constituent name
    real(r8)    , intent(in ) :: latvals(:) ! lat in degrees (ncol)
    real(r8)    , intent(in ) :: lonvals(:) ! lon in degrees (ncol)
    logical     , intent(in ) :: mask   (:) ! Only initialize where .true.
    real(r8)    , intent(out) :: q    (:,:) ! kg tracer/kg dry air (gcol, plev

    integer k

#ifdef USE_UNICON

    do k = 1, size(q, 2)
      if (name == 'ORGawk') then
        where(mask)
         q(:,k) = 0.0_r8
        end where
      else if (name == 'ORGthl') then
       where(mask)
          q(:,k) = 100.0_r8
        end where
      else if (name == 'ORGqto') then
        where (mask)
          q(:,k) = 100.0_r8
        end where
      else if (name == 'ORGuoo') then
        where (mask)
          q(:,k) = 100.0_r8
        end where
      else if (name == 'ORGvoo') then
        where (mask)
          q(:,k) = 100.0_r8
        end where
      end if
    end do
#endif

  end subroutine unicon_init_cnst

  subroutine unicon_cam_tend(dt, state, cam_in, pbuf, ptend, out)

    real(r8)                 , intent(in)  :: dt         ! Time step [s]
    type(physics_state)      , intent(in)  :: state      ! Physics state variables
    type(cam_in_t)           , intent(in)  :: cam_in     ! import state
    type(physics_buffer_desc), pointer     :: pbuf(:)    ! physics buffer
    type(physics_ptend)      , intent(out) :: ptend      ! parameterization tendencies
    type(unicon_out_t)       , intent(out) :: out        ! parameterization outputs

    real(r8) sten_ori(pcols,pver)         !  Tendency of dry static energy [ J / kg / s ]
    real(r8) qvten_ori(pcols,pver)        !  Tendency of water vapor specific humidity [ kg / kg / s ]
    real(r8) qlten_ori(pcols,pver)        !  Tendency of liquid water mixing ratio [ kg / kg / s ]
    real(r8) qiten_ori(pcols,pver)        !  Tendency of ice mixing ratio [ kg / kg / s ]
    real(r8) trten_ori(pcols,pver,PCNST)  !  Tendency of tracers [ # / kg / s, kg / kg / s ]

    real(r8) slten_pos_inv(pcols,pver)    !
    real(r8) qtten_pos_inv(pcols,pver)    !
    real(r8) uten_pos_inv(pcols,pver)     !
    real(r8) vten_pos_inv(pcols,pver)     !
    real(r8) sten_pos_inv(pcols,pver)     !
    real(r8) qvten_pos_inv(pcols,pver)    !
    real(r8) qlten_pos_inv(pcols,pver)    !
    real(r8) qiten_pos_inv(pcols,pver)    !
    real(r8) trten_pos_inv(pcols,pver,PCNST)

    integer iend
    integer lchnk
    integer itim

    ! fields in physics buffer
    real(r8), pointer, dimension(:) :: & ! (mix)
       sgh30       ! Standard deviation of topography at 30s [m]

    real(r8), pointer, dimension(:,:) :: & ! (mix,mkx)
       ast0_inv    ! Physical stratus fraction [ fraction ]

    real(r8), pointer, dimension(:,:) :: & ! (mix,mkx+1)
       tke0_inv,    &! TKE at the interface [ m2/s2 ]
       bprod0_inv    ! Buoyancy production at the interface [ m2/s3 ]

    integer(i4), pointer, dimension(:) :: & ! (mix)
       kpblh_inv       ! Layer index with PBL top in it or at the base interface

    real(r8), pointer, dimension(:) :: & ! (mix)
       pblh,          &! PBL top height [m]
       went,          &! Entrainment rate at the PBL top interface directly from UW PBL scheme [ m / s ]. went = 0 for STL.
       cush,          &! Cumulus top height [ m ]
       cushavg,       &! Mean cumulus top height weighted by updraft mass flux at surface [ m ]
       cuorg,         &! Convective organization parameter [ 0-1 ]

       awk_PBL,       &! Wake area within PBL [ 0 - 1 ]
       delta_thl_PBL, &! Diff of thl between off-wake region and grid-mean value averaged over the PBL [ K ]
       delta_qt_PBL,  &! Diff of qt  between off-wake region and grid-mean value averaged over the PBL [ kg/kg ]
       delta_u_PBL,   &! Diff of u   between off-wake region and grid-mean value averaged over the PBL [ m/s ]
       delta_v_PBL     ! Diff of v   between off-wake region and grid-mean value averaged over the PBL [ m/s ]

    real(r8), pointer, dimension(:,:) :: & ! (mix,pcnst)
       delta_tr_PBL ! Diff of tr  between off-wake region and grid-mean value avg over the PBL [ kg/kg, #/kg ]

    real(r8), dimension(pcols,pver) :: & ! (mix,mkx)
       cu_cmfum  ! The mass involved in the updraft buoyancy sorting at the previous time step [ kg/s/m2 ]

    real(r8), pointer, dimension(:,:) :: & ! (mix,mkx)
       cu_cmfr, &! The detrained mass from convective up and downdraft at the previous time step [ kg/s/m2 ]
       cu_thlr, &! Mass-flux wghted mean 'thl' of detrained mass from conv up and downdraft at prev step [ K ]
       cu_qtr,  &! Mass-flux wghted mean 'qt'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]
       cu_ur,   &! Mass-flux wghted mean 'u'   of detrained mass from conv up and downdraft at prev step [ m/s ]
       cu_vr,   &! Mass-flux wghted mean 'v'   of detrained mass from conv up and downdraft at prev step [ m/s ]
       cu_qlr,  &! Mass-flux wghted mean 'ql'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]
       cu_qir,  &! Mass-flux wghted mean 'qi'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]
       cmfr_det,&! The detrained mass from convective up and downdraft at the previous time step [ kg/s/m2 ]
       qlr_det, &! Mass-flux wghted mean 'ql'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]
       qir_det   ! Mass-flux wghted mean 'qi'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]

    real(r8), pointer, dimension(:,:,:) :: & ! (mix,mkx,pcnst)
       cu_trr    ! Mass-flux wghted mean 'tr'  of detrained mass from conv up and downdraft at prev step [ kg/kg ]

    real(r8), dimension(pcols,pver) :: & ! (mix,mkx)
       cu_cmfrd,&! The amount of detrained mass from convective downdraft at the previous time step [ kg/s/m2 ]
       cu_thlrd,&! Mass-flux wghted mean 'thl' of detrained mass from conv downdraft at prev step [ K ]
       cu_qtrd, &! Mass-flux wghted mean 'qt'  of detrained mass from conv downdraft at prev step [ kg/kg ]
       cu_urd,  &! Mass-flux wghted mean 'u'   of detrained mass from conv downdraft at prev step [ m/s ]
       cu_vrd,  &! Mass-flux wghted mean 'v'   of detrained mass from conv downdraft at prev step [ m/s ]
       cu_qlrd, &! Mass-flux wghted mean 'ql'  of detrained mass from conv downdraft at prev step [ kg/kg ]
       cu_qird   ! Mass-flux wghted mean 'qi'  of detrained mass from conv downdraft at prev step [ kg/kg ]

    real(r8), dimension(pcols,pver,PCNST) :: & ! (mix,mkx,pcnst)
       cu_trrd   ! Mass-flux wghted mean 'tr'  of detrained mass from conv downdraft at prev step [ kg/kg ]

    real(r8), pointer, dimension(:,:) :: & ! (mix,mkx)
       shfrc,     &! Convective updraft fractional area
       icwmrsh,   &! In-cloud LWC+IWC within convective updraft
       rprdsh,    &! Prod of rain+snow by lateral expels of cumulus condensate [ kg / kg / s ]
       evapc_inv, &! Evaporation rate of convective precipitation within environment [ kg/kg/s ]
       am_evp_st_inv,   &!  Evaporation area of stratiform precipitation [fraction]
       evprain_st_inv,  &!  Grid-mean evaporation rate of stratiform rain [kg/kg/s] >= 0.
       evpsnow_st_inv    !  Grid-mean evaporation rate of stratiform snow [kg/kg/s] >= 0.

    real(r8), pointer, dimension(:) :: & ! (mix)
       precip,    &!  Precipitation flux at surface in flux unit [ m / s ]
       snow        !  Snow flux at surface in flux unit [ m / s ]

    real(r8) ps0       (pcols,0:pver)       ! Environmental pressure at full sigma levels
    real(r8) zs0       (pcols,0:pver)       ! Environmental height at full sigma levels
    real(r8) p0        (pcols,  pver)       ! Environmental pressure at half sigma levels
    real(r8) z0        (pcols,  pver)       ! Environmental height at half sigma levels
    real(r8) dp0       (pcols,  pver)       ! Environmental layer pressure thickness
    real(r8) dpdry0    (pcols,  pver)       ! Environmental layer dry pressure thickness
    real(r8) u0        (pcols,  pver)       ! Environmental zonal wind
    real(r8) v0        (pcols,  pver)       ! Environmental meridional wind
    real(r8) qv0       (pcols,  pver)       ! Environmental specific humidity
    real(r8) ql0       (pcols,  pver)       ! Environmental liquid water mixing ratio
    real(r8) qi0       (pcols,  pver)       ! Environmental ice mixing ratio
    real(r8) tr0       (pcols,  pver,PCNST) ! Environmental tracers [ #/kg, kg/kg ]
    real(r8) t0        (pcols,  pver)       ! Environmental temperature
    real(r8) s0        (pcols,  pver)       ! Environmental dry static energy
    real(r8) ast0      (pcols,  pver)       ! Physical stratiform cloud fraction [ fraction ]
    real(r8) tke0      (pcols,0:pver)       ! TKE [ m2/s2 ]
    real(r8) bprod0    (pcols,0:pver)       ! Buoyancy production [ m2/s3 ]
    real(r8) am_evp_st (pcols,  pver)       ! Evaporation area of stratiform precipitation [fraction]
    real(r8) evprain_st(pcols,  pver)       ! Grid-mean evaporation rate of stratiform rain [kg/kg/s] >= 0.
    real(r8) evpsnow_st(pcols,  pver)       ! Grid-mean evaporation rate of stratiform snow [kg/kg/s] >= 0.

    integer(i4) kpblh(pcols)               !  Layer index with PBL top in it or at the base interface

    real(r8) am_u(pcols,pver)               !  Convective updraft fractional area
    real(r8) qlm_u(pcols,pver)              !  In-cloud LWC within convective updraft [ kg / kg ]
    real(r8) qim_u(pcols,pver)              !  In-cloud IWC within convective updraft [ kg / kg ]

    real(r8) am_d(pcols,pver)               !  Convective downdraft fractional area
    real(r8) qlm_d(pcols,pver)              !  In-cloud LWC within downdraft updraft [ kg / kg ]
    real(r8) qim_d(pcols,pver)              !  In-cloud IWC within downdraft updraft [ kg / kg ]

    real(r8) cmf_u(pcols,0:pver)            !  Upward     convective mass flux at the interface [ kg / s / m2 ]
    real(r8) slflx(pcols,0:pver)            !  Net upward convective flux of liquid static energy [ J / s / m2 ]
    real(r8) qtflx(pcols,0:pver)            !  Net upward convective flux of total specific humidity [ kg / s / m2 ]

    real(r8) qvten(pcols,pver)              !  Tendency of water vapor specific humidity [ kg / kg / s ]
    real(r8) qlten(pcols,pver)              !  Tendency of liquid water mixing ratio [ kg / kg / s ]
    real(r8) qiten(pcols,pver)              !  Tendency of ice mixing ratio [ kg / kg / s ]
    real(r8) trten(pcols,pver,PCNST)        !  Tendency of tracers [ # / kg / s, kg / kg / s ]

    real(r8) sten(pcols,pver)               !  Tendency of dry static energy [ J / kg / s ]
    real(r8) uten(pcols,pver)               !  Tendency of zonal wind [ m / s / s ]
    real(r8) vten(pcols,pver)               !  Tendency of meridional wind [ m / s / s ]

    real(r8) qrten(pcols,pver)              !  Production rate of rain by lateral expels of cumulus condensate [ kg / kg / s ]
    real(r8) qsten(pcols,pver)              !  Production rate of snow by lateral expels of cumulus condensate [ kg / kg / s ]

    real(r8) evapc(pcols,pver)              !  Evaporation rate of convective precipitation within environment [ kg/kg/s ]

    real(r8) rqc  (pcols,pver)              !  Production rate of detrained LWC+IWC [ kg / kg / s ] > 0
    real(r8) rqc_l(pcols,pver)              !  Production rate of detrained LWC     [ kg / kg / s ] > 0
    real(r8) rqc_i(pcols,pver)              !  Production rate of detrained IWC     [ kg / kg / s ] > 0
    real(r8) rnc_l(pcols,pver)              !  Production rate of detrained droplet number of cloud liquid droplets [ # / kg / s ] > 0
    real(r8) rnc_i(pcols,pver)              !  Production rate of detrained droplet number of cloud    ice droplets [ # / kg / s ] > 0

    real(r8) cnt(pcols)                    !  Cloud top  interface index ( ki = kpen )
    real(r8) cnb(pcols)                    !  Cloud base interface index ( ki = krel-1 )

    real(r8) pdel0  (pcols,pver)            !  Environmental pressure thickness ( either dry or moist ) [ Pa ]
    real(r8) trmin                       !  Minimum concentration of tracer [ # / kg ]

    real(r8) cmf_det(pcols,pver)            ! Det mass flux from convective updraft (not from environment) and downdraft [kg/s/m^2]
    real(r8) ql_det (pcols,pver)            ! Det conv LWC from convective updraft (not from environment) and downdraft [kg/s/m^2]
    real(r8) qi_det (pcols,pver)            ! Det conv IWC from convective updraft (not from environment) and downdraft [kg/s/m^2]

    ! For prognostically updated state variables

    real(r8) qv0_c(pcols,pver)               !  Environmental specific humidity
    real(r8) ql0_c(pcols,pver)               !  Environmental liquid water mixing ratio
    real(r8) qi0_c(pcols,pver)               !  Environmental ice mixing ratio
    real(r8) t0_c (pcols,pver)               !  Environmental temperature
    real(r8) s0_c (pcols,pver)               !  Environmental dry static energy
    real(r8) tr0_c(pcols,pver,PCNST)         !  Environmental tracers [ # / kg, kg / kg ]

    ! Layer index variables
    integer k                             !  Vertical index for local fields
    integer k_inv                         !  Vertical index for incoming fields
    integer mt                            !  Tracer index [ no ]
    integer m

    ! For aerosol tendency output
    character(30) varname

    logical lq(PCNST)

    ! --------- !
    ! Main body !
    ! --------- !

#ifdef USE_UNICON
    iend  = state%ncol
    lchnk = state%lchnk

    ! Associate pointers with physics buffer fields
    itim = pbuf_old_tim_idx()
    call pbuf_get_field(pbuf, sgh30_idx        , sgh30         )
    call pbuf_get_field(pbuf, ast_idx          , ast0_inv      , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, tke_idx          , tke0_inv      )
    call pbuf_get_field(pbuf, bprod_idx        , bprod0_inv    )
    call pbuf_get_field(pbuf, kpblh_idx        , kpblh_inv     )
    call pbuf_get_field(pbuf, pblh_idx         , pblh          )
    call pbuf_get_field(pbuf, went_idx         , went          )
    call pbuf_get_field(pbuf, cush_idx         , cush          , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, cushavg_idx      , cushavg       , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, cuorg_idx        , cuorg         , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, awk_PBL_idx      , awk_PBL       , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, delta_thl_PBL_idx, delta_thl_PBL , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, delta_qt_PBL_idx , delta_qt_PBL  , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, delta_u_PBL_idx  , delta_u_PBL   , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, delta_v_PBL_idx  , delta_v_PBL   , start=[1,    itim], kount=[pcols,           1])
    call pbuf_get_field(pbuf, delta_tr_PBL_idx , delta_tr_PBL  , start=[1,1,  itim], kount=[pcols,     pcnst,1])
    call pbuf_get_field(pbuf, cu_cmfr_idx      , cu_cmfr       , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_thlr_idx      , cu_thlr       , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_qtr_idx       , cu_qtr        , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_ur_idx        , cu_ur         , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_vr_idx        , cu_vr         , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_qlr_idx       , cu_qlr        , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_qir_idx       , cu_qir        , start=[1,1,  itim], kount=[pcols,pver,      1])
    call pbuf_get_field(pbuf, cu_trr_idx       , cu_trr        , start=[1,1,1,itim], kount=[pcols,pver,pcnst,1])
    call pbuf_get_field(pbuf, shfrc_idx        , shfrc         )
    call pbuf_get_field(pbuf, icwmrsh_idx      , icwmrsh       )
    call pbuf_get_field(pbuf, rprdsh_idx       , rprdsh        )
    call pbuf_get_field(pbuf, prec_sh_idx      , precip        )
    call pbuf_get_field(pbuf, snow_sh_idx      , snow          )
    call pbuf_get_field(pbuf, nevapr_shcu_idx  , evapc_inv     )
    call pbuf_get_field(pbuf, am_evp_st_idx    , am_evp_st_inv )
    call pbuf_get_field(pbuf, evprain_st_idx   , evprain_st_inv)
    call pbuf_get_field(pbuf, evpsnow_st_idx   , evpsnow_st_inv)
    call pbuf_get_field(pbuf, cmfr_det_idx     , cmfr_det      )
    call pbuf_get_field(pbuf, qlr_det_idx      , qlr_det       )
    call pbuf_get_field(pbuf, qir_det_idx      , qir_det       )

    ! Reverse variables defined at the layer mid-point
    do k = 1, pver
      k_inv               = pver + 1 - k
      p0        (:iend,k) = state%pmid    (:iend,k_inv)
      u0        (:iend,k) = state%u       (:iend,k_inv)
      v0        (:iend,k) = state%v       (:iend,k_inv)
      z0        (:iend,k) = state%zm      (:iend,k_inv)
      dp0       (:iend,k) = state%pdel    (:iend,k_inv)
      dpdry0    (:iend,k) = state%pdeldry (:iend,k_inv)
      qv0       (:iend,k) = state%q       (:iend,k_inv,1)
      ql0       (:iend,k) = state%q       (:iend,k_inv,ixcldliq)
      qi0       (:iend,k) = state%q       (:iend,k_inv,ixcldice)
      t0        (:iend,k) = state%t       (:iend,k_inv)
      s0        (:iend,k) = state%s       (:iend,k_inv)
      ast0      (:iend,k) = ast0_inv      (:iend,k_inv)
      am_evp_st (:iend,k) = am_evp_st_inv (:iend,k_inv)
      evprain_st(:iend,k) = evprain_st_inv(:iend,k_inv)
      evpsnow_st(:iend,k) = evpsnow_st_inv(:iend,k_inv)
      do mt = 1, pcnst
        tr0(:iend,k,mt)  = state%q(:iend,k_inv,mt)
      end do
    end do

   ! Reverse variables defined at the interfaces
    do k = 0, pver
      k_inv           = pver + 1 - k
      ps0   (:iend,k) = state%pint(:iend,k_inv)
      zs0   (:iend,k) = state%zi(:iend,k_inv)
      tke0  (:iend,k) = tke0_inv(:iend,k_inv)
      bprod0(:iend,k) = bprod0_inv(:iend,k_inv)
    end do

    ! Reverse the index of ambiguous layer
    kpblh(:iend) = pver + 1 - kpblh_inv(:iend)

    call compute_unicon(pcols       , pver        , iend      , pcnst    , dt      ,           &
                        ps0       , zs0        , p0        , z0       , dp0     , dpdry0  , &
                        t0        , qv0        , ql0       , qi0      , tr0     ,           &
                        u0        , v0         , ast0      , tke0     , bprod0  ,           &
                        kpblh     , pblh       , went      ,                                &
                        cam_in%cflx(:,1), cam_in%shf, cam_in%wsx, cam_in%wsy, cam_in%cflx , &
                        cam_in%landfrac, sgh30   ,                                          &
                        am_evp_st , evprain_st , evpsnow_st,                                &
                        cush      , cushavg    , cuorg     ,                                &
                        awk_PBL                , delta_thl_PBL        , delta_qt_PBL      , &
                        delta_u_PBL            , delta_v_PBL          , delta_tr_PBL      , &
                        cu_cmfum  , cu_cmfr    , cu_thlr   , cu_qtr   , cu_ur   , cu_vr   , &
                        cu_qlr    , cu_qir     , cu_trr    ,                                &
                        cu_cmfrd  , cu_thlrd   , cu_qtrd   , cu_urd   , cu_vrd  ,           &
                        cu_qlrd   , cu_qird    , cu_trrd   ,                                &
                        am_u      , qlm_u      , qim_u     ,                                &
                        am_d      , qlm_d      , qim_d     ,                                &
                        cmf_u     , slflx      , qtflx     ,                                &
                        qvten     , qlten      , qiten     , trten    ,                     &
                        sten      , uten       , vten      ,                                &
                        qrten     , qsten      ,                                            &
                        rqc_l     , rqc_i      , rqc       , rnc_l    , rnc_i   ,           &
                        out%rliq  , precip     , snow      , evapc    ,                     &
                        cnt       , cnb        , cmf_det   , ql_det   , qi_det  ,           &
                        lchnk )

    ! Initialize output ptend
    lq(:) = .true.
    call physics_ptend_init(ptend, state%psetcols, 'unicon', ls=.true., lu=.true., lv=.true., lq=lq)

    ! ----------------------------------------------------- !
    ! Treatment of reserved liquid-ice water for CAM        !
    ! All the reserved condensate is converted into liquid. !
    ! Jan.04.2012. Relocated from below to here to prevent  !
    !              energy and water conservation error.     !
    !              Also add corresponding tendency of cloud !
    !              droplet number concentration.            !
    ! Note that since 'positive_moisture' will convert      !
    ! negative 'ql,qi,nl,ni' into positive, below may cause !
    ! larger 'ql,qi' and smaller 'qv' than it should be if  !
    ! 'ql0,qi0' further below becomes negative.             !
    ! This feature is inevitable at the current stage.      !
    ! Note that in future, I can impose consistency between !
    ! positive_moisture and positive_tracer for nl,ni.      !
    ! ----------------------------------------------------- !

    qlten(:iend,:pver) = qlten(:iend,:pver) - rqc_l(:iend,:pver)
    qiten(:iend,:pver) = qiten(:iend,:pver) - rqc_i(:iend,:pver)
    sten (:iend,:pver) =  sten(:iend,:pver) - (xls - xlv) * rqc_i(:iend,:pver)
    trten(:iend,:pver,ixnumliq) = trten(:iend,:pver,ixnumliq) - rnc_l(:iend,:pver)
    trten(:iend,:pver,ixnumice) = trten(:iend,:pver,ixnumice) - rnc_i(:iend,:pver)

    ! --------------------------------------------- !
    ! Prevent negative cloud condensate and tracers !
    ! --------------------------------------------- !

    qv0_c(:iend,:pver) = qv0(:iend,:pver) + qvten(:iend,:pver) * dt
    ql0_c(:iend,:pver) = ql0(:iend,:pver) + qlten(:iend,:pver) * dt
    qi0_c(:iend,:pver) = qi0(:iend,:pver) + qiten(:iend,:pver) * dt
     t0_c(:iend,:pver) =  t0(:iend,:pver) +  sten(:iend,:pver) * dt / cp
     s0_c(:iend,:pver) =  s0(:iend,:pver) +  sten(:iend,:pver) * dt
    do mt = 1, pcnst
      tr0_c(:iend,:pver,mt) = tr0(:iend,:pver,mt) + trten(:iend,:pver,mt)*dt
    end do

    ! Note : Since 'positive_moisture' will only perform condensation not the evaporation,
    !        we don't need to modify corresponding 'nl,ni'.
    !        Thus, current version is completely OK.

    sten_ori (:iend,:pver) =  sten(:iend,:pver)
    qvten_ori(:iend,:pver) = qvten(:iend,:pver)
    qlten_ori(:iend,:pver) = qlten(:iend,:pver)
    qiten_ori(:iend,:pver) = qiten(:iend,:pver)
    do mt = 1, pcnst
      trten_ori(:iend,:pver,mt) = trten(:iend,:pver,mt)
    end do

    call positive_moisture( &
      cp, xlv, xls, pcols, iend, &
      pver, dt, qmin(1), qmin(ixcldliq), qmin(ixcldice), &
      dp0, qv0_c, ql0_c, qi0_c, t0_c, &
      s0_c, qvten, qlten, qiten, sten)

    do mt = 1, pcnst
      if (cnst_get_type_byind(mt) == 'wet') then
        pdel0(:iend,:pver) = dp0(:iend,:pver)
      else
        pdel0(:iend,:pver) = dpdry0(:iend,:pver)
      endif

      if (cnst_is_mam_num(mt)) then
        trmin = 1.e-5_r8
      else
        trmin = qmin(mt)
      end if
      call positive_tracer(pcols, iend, pver, dt, trmin, pdel0, tr0_c(:,:,mt), trten(:,:,mt))
    end do

    ! ----------------------------------------------------- !
    ! Treatment of reserved liquid-ice water for CAM        !
    ! All the reserved condensate is converted into liquid. !
    ! ----------------------------------------------------- !

    ! Important : While below treatment looks different from what is being done in the uwshcu.F90,
    !             below is actually the same as what is being done in the uwshcu.F90. The process
    !             used below is (1) convert detrained ice into liquid, which involves melting cooling,
    !             and then (2) subtract this total detrained liquid + ice from the grid-mean qlten.
    ! Sep.08.2010. In the new scheme, it will be 'rqc_l = rqc_i = 0'. Thus, below block does nothing.
    !              But it is no harm to keep below block.
    ! Jan.04.2012. Below block will be used again for many reasons : (1) inctease TGCLDLWP, (2) reduce PREH20,
    !              (2) reduce too strong SWCF over the far eastern equatorial Pacific.
    !              Currently, below produces conservation error of both energy and water.
    ! Jan.04.2012. Conservation errors of energy and moisture dissappears if I locate below block
    !              before applying 'positive_moisture'. Thus, I relocated below block above.

    ! qlten(:iend,:mkx) = qlten(:iend,:mkx) + rqc_i(:iend,:mkx)
    ! qiten(:iend,:mkx) = qiten(:iend,:mkx) - rqc_i(:iend,:mkx)
    ! sten(:iend,:mkx)  =  sten(:iend,:mkx) - ( xls - xlv ) * rqc_i(:iend,:mkx)
    ! qlten(:iend,:mkx) = qlten(:iend,:mkx) - rqc_l(:iend,:mkx) - rqc_i(:iend,:mkx)
    ! trten(:iend,:mkx,ixnumliq) = trten(:iend,:mkx,ixnumliq) - rnc_l(:iend,:mkx)
    ! trten(:iend,:mkx,ixnumice) = trten(:iend,:mkx,ixnumice) - rnc_i(:iend,:mkx)

    ! --------------------------- !
    ! Reverse vertical coordinate !
    ! --------------------------- !

    ! Reverse cloud top/base interface indices
    out%cnt(:iend) = real(pver,r8) + 1 - cnt(:iend)
    out%cnb(:iend) = real(pver,r8) + 1 - cnb(:iend)

    ! Reverse variables defined at the interfaces
    do k = 0, pver
      k_inv                  = pver + 1 - k
      out%cmfmc(:iend,k_inv) = cmf_u(:iend,k) !  Updraft mass flux at top of layer
      out%slflx(:iend,k_inv) = slflx(:iend,k) !  Net liquid static energy flux
      out%qtflx(:iend,k_inv) = qtflx(:iend,k) !  Net total water flux
    end do

    ! Reverse variables defined at the layer mid-point
    do k = 1, pver
      k_inv = pver + 1 - k
      ptend%q(:iend,k_inv,1)        = qvten(:iend,k) ! Convective tendency of specific humidity
      ptend%q(:iend,k_inv,ixcldliq) = qlten(:iend,k) ! Convective tendency of liquid water mixing ratio
      ptend%q(:iend,k_inv,ixcldice) = qiten(:iend,k) ! Convective tendency of ice mixing ratio
      ptend%s(:iend,k_inv)          = sten (:iend,k) ! Convective tendency of static energy
      ptend%u(:iend,k_inv)          = uten (:iend,k) ! Convective tendency of zonal wind
      ptend%v(:iend,k_inv)          = vten (:iend,k) ! Convective tendency of meridional wind
      out%rqc(:iend,k_inv)          = rqc  (:iend,k) ! Convective tendency of reserved (suspended) liquid + ice water

      ! These quantities are being put into physics buffer fields that were meant for shallow convection.
      ! This should be fixed.
      shfrc    (:iend,k_inv) = am_u (:iend,k)
      icwmrsh  (:iend,k_inv) = qlm_u(:iend,k) + qim_u(:iend,k)
      rprdsh   (:iend,k_inv) = qrten(:iend,k) + qsten(:iend,k)
      evapc_inv(:iend,k_inv) = evapc(:iend,k)       ! Evaporation rate of convective precipitation within environment

      do mt = 2, pcnst
        if (mt /= ixcldliq .and. mt /= ixcldice) then
          ptend%q(:iend,k_inv,mt) = trten(:iend,k,mt)
        end if
      end do

      ! Additional diagnostic output associated with positive tracer
      sten_pos_inv (:iend,k_inv) = sten (:iend,k) -  sten_ori(:iend,k)
      qvten_pos_inv(:iend,k_inv) = qvten(:iend,k) - qvten_ori(:iend,k)
      qlten_pos_inv(:iend,k_inv) = qlten(:iend,k) - qlten_ori(:iend,k)
      qiten_pos_inv(:iend,k_inv) = qiten(:iend,k) - qiten_ori(:iend,k)
      qtten_pos_inv(:iend,k_inv) = qvten_pos_inv(:iend,k_inv) + qlten_pos_inv(:iend,k_inv) &
                                                              + qiten_pos_inv(:iend,k_inv)
      slten_pos_inv(:iend,k_inv) = sten_pos_inv(:iend,k_inv) - xlv * qlten_pos_inv(:iend,k_inv) &
                                                             - xls * qiten_pos_inv(:iend,k_inv)
      uten_pos_inv(:iend,k_inv)  = 0._r8
      vten_pos_inv(:iend,k_inv)  = 0._r8
      do mt = 1, pcnst
        trten_pos_inv(:iend,k_inv,mt)  = trten(:iend,k,mt) - trten_ori(:iend,k,mt)
      end do

      ! Save detrainment variables to pbuf for use in the cloud macrophysics
      cmfr_det(:iend,k_inv) = cmf_det(:iend,k)
      qlr_det (:iend,k_inv) =  ql_det(:iend,k)
      qir_det (:iend,k_inv) =  qi_det(:iend,k)
    end do

    call outfld('slten_pos_SP' , slten_pos_inv, pcols, lchnk)
    call outfld('qtten_pos_SP' , qtten_pos_inv, pcols, lchnk)
    call outfld('uten_pos_SP'  ,  uten_pos_inv, pcols, lchnk)
    call outfld('vten_pos_SP'  ,  vten_pos_inv, pcols, lchnk)
    call outfld('sten_pos_SP'  ,  sten_pos_inv, pcols, lchnk)
    call outfld('qvten_pos_SP' , qvten_pos_inv, pcols, lchnk)
    call outfld('qlten_pos_SP' , qlten_pos_inv, pcols, lchnk)
    call outfld('qiten_pos_SP' , qiten_pos_inv, pcols, lchnk)
    call outfld('nlten_pos_SP' , trten_pos_inv(:,:,ixnumliq), pcols, lchnk)
    call outfld('niten_pos_SP' , trten_pos_inv(:,:,ixnumice), pcols, lchnk)
    call outfld('CMFR_DET'     , cmfr_det,      pcols, lchnk)
    call outfld('QLR_DET'      , qlr_det,       pcols, lchnk)
    call outfld('QIR_DET'      , qir_det,       pcols, lchnk)

    do m = 1, pcnst
      if (cnst_is_mam_num(m) .or. cnst_is_mam_mmr(m)) then
        varname = trim(cnst_name(m)) // '_pos_SP'
        call outfld(trim(varname), trten_pos_inv(:,:,m), pcols, lchnk)
      end if
    end do
#endif

  end subroutine unicon_cam_tend

  subroutine unicon_cam_org_diags(state, pbuf)

    ! ------------------------------------------------------------------------
    ! Insert the organization-related heterogeneities computed inside the
    ! UNICON into the tracer arrays here before performing advection.
    ! This is necessary to prevent any modifications of organization-related
    ! heterogeneities by non convection-advection process, such as
    ! dry and wet deposition of aerosols, MAM, etc.
    ! Again, note that only UNICON and advection schemes are allowed to
    ! changes to organization at this stage, although we can include the
    ! effects of other physical processes in future.
    ! ------------------------------------------------------------------------

    type(physics_state), intent(inout) :: state
    type(physics_buffer_desc), pointer :: pbuf(:)

    real(r8), pointer, dimension(:) :: awk_PBL
    real(r8), pointer, dimension(:) :: delta_thl_PBL
    real(r8), pointer, dimension(:) :: delta_qt_PBL
    real(r8), pointer, dimension(:) :: delta_u_PBL
    real(r8), pointer, dimension(:) :: delta_v_PBL

    integer i, itim, ncol

#ifdef USE_UNICON
    ncol = state%ncol

    itim = pbuf_old_tim_idx()
    call pbuf_get_field(pbuf, awk_PBL_idx      , awk_PBL      , start=[1,itim], kount=[pcols,1])
    call pbuf_get_field(pbuf, delta_thl_PBL_idx, delta_thl_PBL, start=[1,itim], kount=[pcols,1])
    call pbuf_get_field(pbuf, delta_qt_PBL_idx , delta_qt_PBL , start=[1,itim], kount=[pcols,1])
    call pbuf_get_field(pbuf, delta_u_PBL_idx  , delta_u_PBL  , start=[1,itim], kount=[pcols,1])
    call pbuf_get_field(pbuf, delta_v_PBL_idx  , delta_v_PBL  , start=[1,itim], kount=[pcols,1])

    do i = 1, ncol
      state%q(i,:,awk_cnst_ind) = awk_PBL(i)
      ! Add a constant offset of '100._r8' to 'delta_xxx' variables to be
      ! consistent with the reading sentence of unicon.F90 and so to prevent
      ! negative tracers.
      state%q(i,:,thl_cnst_ind) = max(0.0_r8, delta_thl_PBL(i) + 100.0_r8)
      state%q(i,:, qt_cnst_ind) = max(0.0_r8,  delta_qt_PBL(i) + 100.0_r8)
      state%q(i,:,  u_cnst_ind) = max(0.0_r8,   delta_u_PBL(i) + 100.0_r8)
      state%q(i,:,  v_cnst_ind) = max(0.0_r8,   delta_v_PBL(i) + 100.0_r8)
    end do
#endif

  end subroutine unicon_cam_org_diags

end module unicon_cam
